Basic statistics for health research - week 10, ex. 1

Prep and setup

library(nlme)
library(lattice)
library(lindia)
library(mgcv)
library(tidyverse)
library(descr)
library(patchwork)

theme_set(theme_minimal())

pretty_estimates <- function(m) { # m: model object
    cbind(Est. = m$coefficients, confint(m))
}

Load data

lang <- read_delim("http://publicifsv.sund.ku.dk/~lts/basal/data/varicer.txt", delim = " ") %>% 
    mutate(patient = factor(patient),
           treat = factor(treat, levels = c("P", "A")))

lang5 <- lang %>% 
    filter(minutter <= 5)

Q1

ggplot(lang, aes(x = minutter, y = pressure, colour = patient)) +
    geom_line(alpha = 0.5, show.legend = FALSE) +
    geom_point(size = 1, show.legend = FALSE) +
    facet_wrap(~ treat, nrow = 1) 

lang %>% 
    group_by(treat) %>% 
    summarise(n_patients = n_distinct(patient))

## # A tibble: 2 x 2
##   treat n_patients
##   <fct>      <int>
## 1 P              9
## 2 A              8

lang %>% 
    group_by(patient) %>% 
    summarise(n_obs = n())

## # A tibble: 17 x 2
##    patient n_obs
##    <fct>   <int>
##  1 1          13
##  2 2           9
##  3 3           9
##  4 4          12
##  5 5           9
##  6 6           9
##  7 7           9
##  8 8          11
##  9 9          10
## 10 10         17
## 11 11         17
## 12 12          9
## 13 13         10
## 14 14         14
## 15 15         11
## 16 16         10
## 17 17         10

# Fewer and fewer patients, mainly in placebo arm (not a good sign)
lang %>% 
    group_by(minutter, treat) %>% 
    summarise(n_patients = n_distinct(patient)) %>% 
    ggplot(aes(x = minutter, y = n_patients, colour = treat)) +
        geom_line() +
        coord_cartesian(ylim = c(0, NA))

# Mean VP develops
lang %>% 
    group_by(treat, minutter) %>% 
    summarise(mean_vp = mean(pressure)) %>% 
    ggplot(aes(x = minutter, y = mean_vp, colour = treat)) +
        geom_line() +
        coord_cartesian(ylim = c(0, NA))

Q2

ggplot(lang, aes(x = minutter, y = pressure, colour = patient)) +
    geom_line(alpha = 0.5, show.legend = FALSE) +
    geom_point(size = 1, show.legend = FALSE) +
    facet_wrap(~ treat, nrow = 1) +
    # geom_label(aes(label = patient), ~ filter(., minutter == 0), size = 8/ggplot2::.pt) +
    # facet_wrap(~ patient) +
    coord_cartesian(xlim = c(0, 5)) # Q2

One patient is quite different from the rest in the active arm (suddent drop around 2 minutes in). To identify the person we can add labels to the lines (or plot them individually, see the two commented-out lines above) and we see it patient 13. Let’s remove this patient for now as it seems an outlier of some sort and difficult to model.

Restrict population and time span

lang <- lang %>% 
    filter(patient != 13, minutter <= 5)

Q3

# Plot mean over time by arm
plot_df <- lang %>% 
    group_by(treat, minutter) %>% 
    summarise(mean_vp = mean(pressure),
              sd_vp = sd(pressure), 
              .groups = "drop")

plot_mean_vp <- ggplot(plot_df, aes(x = minutter, y = mean_vp, colour = treat)) +
    geom_line() +
    coord_cartesian(ylim = c(0, NA))

plot_sd_vp <- ggplot(plot_df, aes(x = minutter, y = sd_vp, colour = treat)) +
    geom_line() +
    coord_cartesian(ylim = c(0, NA))

plot_mean_vp /
    plot_sd_vp

This is not expected: values drop too quickly after 3 minutes, due to the low patient count. Otherwise the standard deviation is pretty stable over timer and no too different between groups. This difference it likely that those in the placebo with higher presssures are lost.

At each follow-up time we lose patients, and because those lost patients have higher VP’s, the mean will of course decrease as well.

With respect to randomisation, we need to look at the baseline measurements, where the two groups seem quite comparable. Little differences are to be expected as it’s a little population and randomisation can’t account for random variance. We shouldn’t do a test per se, but if one were to do it it would be a T-test:

t.test(pressure ~ treat, data = filter(lang, minutter == 0))

## 
##  Welch Two Sample t-test
## 
## data:  pressure by treat
## t = -0.60929, df = 12.462, p-value = 0.5533
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -9.687360  5.439741
## sample estimates:
## mean in group P mean in group A 
##        19.23333        21.35714

We could test the difference at a specific time point, say, after 5 minutes:

t.test(pressure ~ treat, data = filter(lang, minutter == 5))

## 
##  Welch Two Sample t-test
## 
## data:  pressure by treat
## t = 0.10015, df = 4.4469, p-value = 0.9246
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -7.013492  7.560158
## sample estimates:
## mean in group P mean in group A 
##        14.53333        14.26000

with(t.test(pressure ~ treat, data = filter(lang, minutter == 5)),
     c(diff(rev(estimate)), conf.int))

## mean in group P                                 
##       0.2733333      -7.0134915       7.5601582

table(filter(lang, minutter == 5)$treat) # very few observations, so non-parametric test perhaps better

## 
## P A 
## 3 5

wilcox.test(pressure ~ treat, data = filter(lang, minutter == 5))

## 
##  Wilcoxon rank sum test
## 
## data:  pressure by treat
## W = 9, p-value = 0.7857
## alternative hypothesis: true location shift is not equal to 0

Anyway, don’t read too much into this as we didn’t account to loss to follow-up. But we do start a table for the results:

result <- tibble(metode = "3d: T-test efter 5 minutter", 
                 behandlingseffekt_5_minutter = "-0.27 (-7.57, 7.01)", 
                 p_value = 0.92)

e)

Two approaches:

• T-test
• ANCOVA model

lang <- lang %>% 
    mutate(dif05 = pressure - baseline)

t.test(dif05 ~ treat, data = filter(lang, minutter == 5)) 

## 
##  Welch Two Sample t-test
## 
## data:  dif05 by treat
## t = 3.8401, df = 4.4813, p-value = 0.01492
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##   2.332268 12.881065
## sample estimates:
## mean in group P mean in group A 
##       0.9666667      -6.6400000

result <- result %>% 
    rbind(c("3e: T-test på ændring 0-5 minutter", "-7.61 (-12.88, -2.33)", 0.015))

Almost nothing changes in placebo group, but rather large change in active arm. We add the signs to results to make them consistent.

ancova <- lm(pressure ~ baseline + treat, data = filter(lang, minutter == 5))
summary(ancova)

## 
## Call:
## lm(formula = pressure ~ baseline + treat, data = filter(lang, 
##     minutter == 5))
## 
## Residuals:
##       1       2       3       4       5       6       7       8 
##  0.2870 -0.0642 -0.2229  1.0531  4.1460 -5.4169 -1.6232  1.8411 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   5.4249     3.6829   1.473   0.2007  
## baseline      0.6714     0.2328   2.884   0.0344 *
## treatA       -5.1968     2.9418  -1.767   0.1376  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.28 on 5 degrees of freedom
## Multiple R-squared:  0.6249, Adjusted R-squared:  0.4748 
## F-statistic: 4.164 on 2 and 5 DF,  p-value: 0.08619

pretty_estimates(ancova)

##                   Est.        2.5 %    97.5 %
## (Intercept)  5.4249496  -4.04223058 14.892130
## baseline     0.6713796   0.07288985  1.269869
## treatA      -5.1967840 -12.75881549  2.365248

result <- result %>% 
    rbind(c("3e: ANCOVA, baseline som kovariat", "-5.20 (-12.76, 2.37)", 0.14))

pred_df <- expand.grid(treat = c("P", "A"),
                       baseline = 10:30) %>% 
    mutate(pressure_pred = predict(ancova, newdata = .))

ggplot(filter(lang, minutter == 5), aes(x = baseline, y = pressure, colour = treat)) +
    geom_point() +
    geom_line(aes(y = pressure_pred), pred_df) +
    labs(x = "Baseline pressure", y = "Pressure at 5 minuttes")

The plot doesn’t convey much information but serves only to show what we model: the treatment effect, which is the (constant!) vertical difference between the curves.

Q4

From now we look at all measurements and not just a specific point in time.

lang <- lang %>% 
    mutate(treatA_minutter = ifelse(treat == "A", minutter, 0),
           minutter_fct = factor(minutter),
           treatA_minutter_fct = factor(treatA_minutter))

xtabs(~ treat + treatA_minutter, data = filter(lang, minutter <= 5))

##      treatA_minutter
## treat  0 0.5  1 1.5  2 2.5  3 3.5  4 4.5  5
##     P 88   0  0   0  0   0  0   0  0   0  0
##     A  7   7  7   7  7   7  7   6  7   6  5

model4 <- lme(pressure ~ minutter_fct + treatA_minutter_fct,
              random = ~ 1 | patient, 
              na.action = na.exclude, 
              data = filter(lang, minutter <= 5))
summary(model4) # disregard the correlation matrix that just clutters everything

## Linear mixed-effects model fit by REML
##  Data: filter(lang, minutter <= 5) 
##        AIC      BIC    logLik
##   699.3604 767.0182 -326.6802
## 
## Random effects:
##  Formula: ~1 | patient
##         (Intercept) Residual
## StdDev:     7.12126  1.63761
## 
## Fixed effects: pressure ~ minutter_fct + treatA_minutter_fct 
##                            Value Std.Error  DF   t-value p-value
## (Intercept)            20.162500 1.8267819 125 11.037169  0.0000
## minutter_fct0.5         0.275554 0.7677241 125  0.358923  0.7203
## minutter_fct1           0.486665 0.7677241 125  0.633906  0.5273
## minutter_fct1.5         1.164443 0.7677241 125  1.516747  0.1319
## minutter_fct2           1.364443 0.7677241 125  1.777257  0.0780
## minutter_fct2.5         1.042221 0.7677241 125  1.357546  0.1771
## minutter_fct3           0.608887 0.7677241 125  0.793107  0.4292
## minutter_fct3.5         0.819998 0.7677241 125  1.068090  0.2875
## minutter_fct4           0.642221 0.7677241 125  0.836525  0.4045
## minutter_fct4.5         1.574489 1.0015045 125  1.572124  0.1184
## minutter_fct5           1.758823 1.1174624 125  1.573944  0.1180
## treatA_minutter_fct0.5 -0.244123 1.1523709 125 -0.211845  0.8326
## treatA_minutter_fct1   -1.012377 1.1523709 125 -0.878517  0.3813
## treatA_minutter_fct1.5 -1.604441 1.1523709 125 -1.392296  0.1663
## treatA_minutter_fct2   -2.761584 1.1523709 125 -2.396437  0.0180
## treatA_minutter_fct2.5 -3.353647 1.1523709 125 -2.910215  0.0043
## treatA_minutter_fct3   -3.906028 1.1523709 125 -3.389558  0.0009
## treatA_minutter_fct3.5 -4.721591 1.1824685 125 -3.992995  0.0001
## treatA_minutter_fct4   -5.339362 1.1523709 125 -4.633371  0.0000
## treatA_minutter_fct4.5 -7.324097 1.3462905 125 -5.440205  0.0000
## treatA_minutter_fct5   -8.423005 1.4685338 125 -5.735656  0.0000
##  Correlation: 
##                        (Intr) mn_0.5 mntt_1 mn_1.5 mntt_2 mn_2.5 mntt_3 mn_3.5
## minutter_fct0.5        -0.120                                                 
## minutter_fct1          -0.120  0.494                                          
## minutter_fct1.5        -0.120  0.494  0.494                                   
## minutter_fct2          -0.120  0.494  0.494  0.494                            
## minutter_fct2.5        -0.120  0.494  0.494  0.494  0.494                     
## minutter_fct3          -0.120  0.494  0.494  0.494  0.494  0.494              
## minutter_fct3.5        -0.120  0.494  0.494  0.494  0.494  0.494  0.494       
## minutter_fct4          -0.120  0.494  0.494  0.494  0.494  0.494  0.494  0.494
## minutter_fct4.5        -0.092  0.379  0.379  0.379  0.379  0.379  0.379  0.379
## minutter_fct5          -0.082  0.340  0.340  0.340  0.340  0.340  0.340  0.340
## treatA_minutter_fct0.5  0.000 -0.657 -0.320 -0.320 -0.320 -0.320 -0.320 -0.320
## treatA_minutter_fct1    0.000 -0.320 -0.657 -0.320 -0.320 -0.320 -0.320 -0.320
## treatA_minutter_fct1.5  0.000 -0.320 -0.320 -0.657 -0.320 -0.320 -0.320 -0.320
## treatA_minutter_fct2    0.000 -0.320 -0.320 -0.320 -0.657 -0.320 -0.320 -0.320
## treatA_minutter_fct2.5  0.000 -0.320 -0.320 -0.320 -0.320 -0.657 -0.320 -0.320
## treatA_minutter_fct3    0.000 -0.320 -0.320 -0.320 -0.320 -0.320 -0.657 -0.320
## treatA_minutter_fct3.5  0.000 -0.312 -0.312 -0.312 -0.312 -0.312 -0.312 -0.640
## treatA_minutter_fct4    0.000 -0.320 -0.320 -0.320 -0.320 -0.320 -0.320 -0.320
## treatA_minutter_fct4.5  0.000 -0.274 -0.274 -0.274 -0.274 -0.274 -0.274 -0.274
## treatA_minutter_fct5    0.000 -0.251 -0.251 -0.251 -0.251 -0.251 -0.251 -0.251
##                        mntt_4 mn_4.5 mntt_5 tA__0. trA__1 tA__1. trA__2 tA__2.
## minutter_fct0.5                                                               
## minutter_fct1                                                                 
## minutter_fct1.5                                                               
## minutter_fct2                                                                 
## minutter_fct2.5                                                               
## minutter_fct3                                                                 
## minutter_fct3.5                                                               
## minutter_fct4                                                                 
## minutter_fct4.5         0.379                                                 
## minutter_fct5           0.340  0.297                                          
## treatA_minutter_fct0.5 -0.320 -0.245 -0.220                                   
## treatA_minutter_fct1   -0.320 -0.245 -0.220  0.487                            
## treatA_minutter_fct1.5 -0.320 -0.245 -0.220  0.487  0.487                     
## treatA_minutter_fct2   -0.320 -0.245 -0.220  0.487  0.487  0.487              
## treatA_minutter_fct2.5 -0.320 -0.245 -0.220  0.487  0.487  0.487  0.487       
## treatA_minutter_fct3   -0.320 -0.245 -0.220  0.487  0.487  0.487  0.487  0.487
## treatA_minutter_fct3.5 -0.312 -0.239 -0.214  0.475  0.475  0.475  0.475  0.475
## treatA_minutter_fct4   -0.657 -0.245 -0.220  0.487  0.487  0.487  0.487  0.487
## treatA_minutter_fct4.5 -0.274 -0.738 -0.215  0.417  0.417  0.417  0.417  0.417
## treatA_minutter_fct5   -0.251 -0.220 -0.756  0.382  0.382  0.382  0.382  0.382
##                        trA__3 tA__3. trA__4 tA__4.
## minutter_fct0.5                                   
## minutter_fct1                                     
## minutter_fct1.5                                   
## minutter_fct2                                     
## minutter_fct2.5                                   
## minutter_fct3                                     
## minutter_fct3.5                                   
## minutter_fct4                                     
## minutter_fct4.5                                   
## minutter_fct5                                     
## treatA_minutter_fct0.5                            
## treatA_minutter_fct1                              
## treatA_minutter_fct1.5                            
## treatA_minutter_fct2                              
## treatA_minutter_fct2.5                            
## treatA_minutter_fct3                              
## treatA_minutter_fct3.5  0.475                     
## treatA_minutter_fct4    0.487  0.475              
## treatA_minutter_fct4.5  0.417  0.406  0.417       
## treatA_minutter_fct5    0.382  0.371  0.382  0.352
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.96134903 -0.51474694  0.01765179  0.47638279  2.40697726 
## 
## Number of Observations: 161
## Number of Groups: 16

Large residuals

resids <- residuals(model4, level = 0) # population level
cov <- extract.lme.cov(model4) 
L <- t(chol(cov)) # transposed cholesky defactorisation of the covariance matrix

plot_df <- tibble(scaled_residuals = solve(L) %*% resids, # solve(L) yields the inversed L
                  fitted_values = fitted(model4,level = 0))

ggplot(plot_df) +
    geom_qq_line(aes(sample = scaled_residuals), colour = "grey60") +
    geom_qq(aes(sample = scaled_residuals), size = 0.5) 

ggplot(plot_df, aes(x = fitted_values, y = scaled_residuals)) +
    geom_jitter(size = 0.5) +
    geom_smooth()

Conditional residuals

qplot(x = fitted(model4,level=1), y = scale(resid(model4,level=1))) +
    geom_hline(yintercept = 0)

qplot(sample = scale(resid(model4)), geom = "qq") +
    geom_qq_line()

qplot(scale(resid(model4)))

b) Quantification of effect

Tungen lige i munden: minutter_fct* correspond to the deviation from the common mean of the placebo group, whereas the treatA_minutter_* corresponds to the difference between the placebo group and active arm. Thus, we can’t give a single estimate but must resort so timepoint-specific estimates because we didn’t make any assumptions about an overall trend from which the two groups could deviate. And we can see significant difference from the 2-minute mark onward.

intervals(model4)

## Approximate 95% confidence intervals
## 
##  Fixed effects:
##                              lower       est.      upper
## (Intercept)             16.5470719 20.1625000 23.7779281
## minutter_fct0.5         -1.2438673  0.2755540  1.7949753
## minutter_fct1           -1.0327561  0.4866651  2.0060864
## minutter_fct1.5         -0.3549784  1.1644429  2.6838642
## minutter_fct2           -0.1549784  1.3644429  2.8838642
## minutter_fct2.5         -0.4772006  1.0422207  2.5616420
## minutter_fct3           -0.9105339  0.6088874  2.1283086
## minutter_fct3.5         -0.6994228  0.8199985  2.3394197
## minutter_fct4           -0.8772006  0.6422207  2.1616420
## minutter_fct4.5         -0.4076123  1.5744894  3.5565911
## minutter_fct5           -0.4527737  1.7588231  3.9704198
## treatA_minutter_fct0.5  -2.5248085 -0.2441235  2.0365615
## treatA_minutter_fct1    -3.2930625 -1.0123774  1.2683076
## treatA_minutter_fct1.5  -3.8851260 -1.6044409  0.6762441
## treatA_minutter_fct2    -5.0422688 -2.7615838 -0.4808988
## treatA_minutter_fct2.5  -5.6343323 -3.3536473 -1.0729623
## treatA_minutter_fct3    -6.1867133 -3.9060282 -1.6253432
## treatA_minutter_fct3.5  -7.0618425 -4.7215906 -2.3813387
## treatA_minutter_fct4    -7.6200466 -5.3393616 -3.0586766
## treatA_minutter_fct4.5  -9.9885732 -7.3240971 -4.6596211
## treatA_minutter_fct5   -11.3294155 -8.4230049 -5.5165942
## attr(,"label")
## [1] "Fixed effects:"
## 
##  Random Effects:
##   Level: patient 
##                   lower    est.   upper
## sd((Intercept)) 4.96722 7.12126 10.2094
## 
##  Within-group standard error:
##    lower     est.    upper 
## 1.446722 1.637610 1.853684

intervals(model4)$fixed[21, ] # estimates at specific time point

##      lower       est.      upper 
## -11.329415  -8.423005  -5.516594

The estimates are much different from the naive ones found earlier (3d) but more like those of 3e, however with much narrower confidence interval. This is because the current analysis included much more data, and this lends information to each specific time point.

We add the results to the table:

result <- result %>% 
    rbind(c("4b: CS-model med minutter som class", "-8.42 (-11.33, -5.53)", 0.0000))

There is a significant difference between the two groups overall.

model4_ml <- lme(pressure ~ minutter_fct + treatA_minutter_fct,
                 random = ~ 1 | patient, 
                 na.action = na.exclude, 
                 data = filter(lang, minutter <= 5),
                 method = "ML")

model4_red <- update(model4_ml, . ~ minutter_fct)

anova(model4_ml, model4_red)

##            Model df      AIC      BIC    logLik   Test L.Ratio p-value
## model4_ml      1 23 723.1673 794.0396 -338.5837                       
## model4_red     2 13 775.6512 815.7095 -374.8256 1 vs 2 72.4839  <.0001

Models with lower AIC and BIC (but larger log-likelihood) are preferred, so modelling the difference between groups improves the model significantly. AIC and BIC account somewhat for the fact that a model with many parameters needs to be a lot better to be preferred.

plot_df <- model4$data %>% 
    mutate(pred4_fix = predict(model4,level=0), # fixed effects
           pred4 = fitted(model4)) # patient-level effects

ggplot(plot_df, aes(x = minutter, y = pred4_fix, colour = treat)) +
    geom_line() +
    geom_point() +
    coord_cartesian(ylim = c(0, NA))

ggplot(plot_df, aes(x = minutter, y = pred4, colour = patient)) +
    geom_line() +
    geom_point() +
    facet_wrap(~ treat, nrow = 1)

Q5

Three things to notice:

• The \(\phi\) (phi) parameter (0.96) is the correlation between two observations in the same patient 1 minute apart. Actual neighbour observations are only 0.5 minutes apart, and they will have a correlation of \(\sqrt{\phi} = \sqrt{0.96} \approx 0.98\).
• Because of the autoregressive nature of the model, there isn’t much room for the random effect of patient (the (Intercept) of the StdDev parameter is close to 0).
• AIC is a lot lower for this model than for the one with an unstructured covariance structure, so we prefer this one and use from now.

model5a <- lme(pressure ~ minutter_fct + treatA_minutter_fct,
               random = ~1 | patient, 
               na.action = na.exclude, 
               correlation = corCAR1(0.9, ~ minutter | patient, FALSE), 
               data = filter(lang, minutter <= 5))
summary(model5a)

## Linear mixed-effects model fit by REML
##  Data: filter(lang, minutter <= 5) 
##        AIC     BIC    logLik
##   611.0626 681.662 -281.5313
## 
## Random effects:
##  Formula: ~1 | patient
##         (Intercept) Residual
## StdDev:  0.01361138 6.854021
## 
## Correlation Structure: Continuous AR(1)
##  Formula: ~minutter | patient 
##  Parameter estimate(s):
##       Phi 
## 0.9639969 
## Fixed effects: pressure ~ minutter_fct + treatA_minutter_fct 
##                            Value Std.Error  DF   t-value p-value
## (Intercept)            20.162500 1.7135087 125 11.766792  0.0000
## minutter_fct0.5         0.305343 0.4346212 125  0.702549  0.4836
## minutter_fct1           0.499881 0.6106497 125  0.818604  0.4146
## minutter_fct1.5         1.161386 0.7430574 125  1.562983  0.1206
## minutter_fct2           1.345410 0.8525002 125  1.578193  0.1170
## minutter_fct2.5         1.007502 0.9470442 125  1.063838  0.2895
## minutter_fct3           0.558767 1.0308598 125  0.542040  0.5888
## minutter_fct3.5         0.754757 1.1064440 125  0.682146  0.4964
## minutter_fct4           0.562132 1.1754375 125  0.478233  0.6333
## minutter_fct4.5         1.724982 1.3304142 125  1.296575  0.1972
## minutter_fct5           1.165922 1.5122688 125  0.770975  0.4422
## treatA_minutter_fct0.5 -0.312211 0.6553980 125 -0.476369  0.6346
## treatA_minutter_fct1   -1.042584 0.9184922 125 -1.135104  0.2585
## treatA_minutter_fct1.5 -1.597455 1.1148100 125 -1.432939  0.1544
## treatA_minutter_fct2   -2.718080 1.2757755 125 -2.130532  0.0351
## treatA_minutter_fct2.5 -3.274290 1.4137012 125 -2.316111  0.0222
## treatA_minutter_fct3   -3.791468 1.5349735 125 -2.470054  0.0149
## treatA_minutter_fct3.5 -4.334556 1.6496563 125 -2.627551  0.0097
## treatA_minutter_fct4   -5.156303 1.7415996 125 -2.960671  0.0037
## treatA_minutter_fct4.5 -7.795249 1.9049122 125 -4.092183  0.0001
## treatA_minutter_fct5   -7.802605 2.0976402 125 -3.719706  0.0003
##  Correlation: 
##                        (Intr) mn_0.5 mntt_1 mn_1.5 mntt_2 mn_2.5 mntt_3 mn_3.5
## minutter_fct0.5        -0.072                                                 
## minutter_fct1          -0.101  0.702                                          
## minutter_fct1.5        -0.123  0.570  0.811                                   
## minutter_fct2          -0.142  0.490  0.698  0.860                            
## minutter_fct2.5        -0.158  0.436  0.620  0.764  0.888                     
## minutter_fct3          -0.173  0.395  0.562  0.693  0.806  0.907              
## minutter_fct3.5        -0.187  0.363  0.517  0.638  0.741  0.834  0.920       
## minutter_fct4          -0.199  0.338  0.481  0.593  0.689  0.775  0.855  0.929
## minutter_fct4.5        -0.196  0.295  0.419  0.517  0.601  0.676  0.745  0.810
## minutter_fct5          -0.190  0.256  0.364  0.449  0.522  0.587  0.647  0.704
## treatA_minutter_fct0.5  0.000 -0.660 -0.461 -0.372 -0.318 -0.281 -0.254 -0.232
## treatA_minutter_fct1    0.000 -0.462 -0.658 -0.531 -0.454 -0.402 -0.362 -0.331
## treatA_minutter_fct1.5  0.000 -0.374 -0.532 -0.656 -0.562 -0.496 -0.448 -0.410
## treatA_minutter_fct2    0.000 -0.321 -0.457 -0.563 -0.655 -0.579 -0.522 -0.477
## treatA_minutter_fct2.5  0.000 -0.284 -0.405 -0.499 -0.580 -0.653 -0.589 -0.539
## treatA_minutter_fct3    0.000 -0.257 -0.366 -0.451 -0.525 -0.591 -0.651 -0.596
## treatA_minutter_fct3.5  0.000 -0.235 -0.334 -0.412 -0.479 -0.540 -0.595 -0.647
## treatA_minutter_fct4    0.000 -0.218 -0.311 -0.383 -0.446 -0.502 -0.553 -0.602
## treatA_minutter_fct4.5  0.000 -0.196 -0.279 -0.344 -0.400 -0.450 -0.497 -0.540
## treatA_minutter_fct5    0.000 -0.175 -0.249 -0.307 -0.357 -0.402 -0.443 -0.482
##                        mntt_4 mn_4.5 mntt_5 tA__0. trA__1 tA__1. trA__2 tA__2.
## minutter_fct0.5                                                               
## minutter_fct1                                                                 
## minutter_fct1.5                                                               
## minutter_fct2                                                                 
## minutter_fct2.5                                                               
## minutter_fct3                                                                 
## minutter_fct3.5                                                               
## minutter_fct4                                                                 
## minutter_fct4.5         0.872                                                 
## minutter_fct5           0.757  0.868                                          
## treatA_minutter_fct0.5 -0.215 -0.186 -0.161                                   
## treatA_minutter_fct1   -0.306 -0.266 -0.229  0.701                            
## treatA_minutter_fct1.5 -0.379 -0.328 -0.284  0.567  0.809                     
## treatA_minutter_fct2   -0.441 -0.383 -0.331  0.486  0.694  0.858              
## treatA_minutter_fct2.5 -0.498 -0.432 -0.373  0.431  0.615  0.760  0.886       
## treatA_minutter_fct3   -0.551 -0.478 -0.413  0.390  0.556  0.687  0.801  0.904
## treatA_minutter_fct3.5 -0.598 -0.519 -0.448  0.356  0.508  0.628  0.732  0.826
## treatA_minutter_fct4   -0.648 -0.562 -0.486  0.331  0.472  0.584  0.681  0.768
## treatA_minutter_fct4.5 -0.582 -0.672 -0.580  0.297  0.424  0.524  0.611  0.690
## treatA_minutter_fct5   -0.519 -0.599 -0.695  0.265  0.378  0.467  0.545  0.615
##                        trA__3 tA__3. trA__4 tA__4.
## minutter_fct0.5                                   
## minutter_fct1                                     
## minutter_fct1.5                                   
## minutter_fct2                                     
## minutter_fct2.5                                   
## minutter_fct3                                     
## minutter_fct3.5                                   
## minutter_fct4                                     
## minutter_fct4.5                                   
## minutter_fct5                                     
## treatA_minutter_fct0.5                            
## treatA_minutter_fct1                              
## treatA_minutter_fct1.5                            
## treatA_minutter_fct2                              
## treatA_minutter_fct2.5                            
## treatA_minutter_fct3                              
## treatA_minutter_fct3.5  0.914                     
## treatA_minutter_fct4    0.850  0.923              
## treatA_minutter_fct4.5  0.763  0.829  0.898       
## treatA_minutter_fct5    0.680  0.739  0.800  0.892
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -1.5330963 -0.9428729 -0.3393333  0.7592998  2.0852057 
## 
## Number of Observations: 161
## Number of Groups: 16

intervals(model5a, which = "fixed")$fixed[21, ] # note we must specify which

##      lower       est.      upper 
## -11.954095  -7.802605  -3.651115

model5a_ml <- update(model5a, method = "ML")
model5a_red <- update(model5a, . ~ minutter_fct, method = "ML")
anova(model5a_ml, model5a_red) # still significantly better model *with* group differences

##             Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## model5a_ml      1 24 620.8233 694.7770 -286.4117                        
## model5a_red     2 14 622.6560 665.7957 -297.3280 1 vs 2 21.83272   0.016

result <- result %>% 
    rbind(c("5: SP(POW)-model med minutter som class", "-7.80 (-12.00, -3.61)", 0.0003))

compute_residuals <- function(model) {
    large_resids <- residuals(model, level = 0)
    cov <- extract.lme.cov(model, lang) 
    L <- t(chol(cov)) 
    
    model$data %>% 
        mutate(large_resids_scaled = as.numeric(solve(L) %*% large_resids),
               conditional_resids_scaled = as.numeric(scale(resid(model))),
               fitted_values = fitted(model, level = 0))
}

ggplot(compute_residuals(model5a), aes(sample = large_resids_scaled)) +
    geom_qq() +
    geom_qq_line()

ggplot(compute_residuals(model5a), aes(x = fitted_values, y= large_resids_scaled)) +
    geom_point() + 
    geom_smooth()

ggplot(compute_residuals(model5a), aes(sample = conditional_resids_scaled)) +
    geom_qq() +
    geom_qq_line()

plot_df <- model5a$data %>% 
    mutate(pred5_fix = predict(model5a, level = 0), # fixed effects
           pred5 = fitted(model5a)) # patient-level effects

ggplot(plot_df, aes(x = minutter, y = pred5_fix, colour = treat)) +
    geom_line() +
    geom_point()

# This plot shows that the predictions for all patients are almost the same unlike model4
ggplot(plot_df, aes(x = minutter, y = pred5, colour = patient)) +
    geom_line() +
    geom_point() +
    facet_wrap(~ treat, nrow = 1)

Q6:

model6a <- lme(pressure ~ minutter + treat:minutter, 
               random = ~ 1 | patient, 
               data = lang,
               correlation = corCAR1(0.9, ~ minutter | patient, FALSE), 
               na.action = na.omit)

summary(model6a)

## Linear mixed-effects model fit by REML
##  Data: lang 
##        AIC      BIC   logLik
##   601.4561 619.8316 -294.728
## 
## Random effects:
##  Formula: ~1 | patient
##         (Intercept) Residual
## StdDev:    4.626822  5.02678
## 
## Correlation Structure: Continuous AR(1)
##  Formula: ~minutter | patient 
##  Parameter estimate(s):
##       Phi 
## 0.9342032 
## Fixed effects: pressure ~ minutter + treat:minutter 
##                     Value Std.Error  DF   t-value p-value
## (Intercept)     20.234116 1.7075505 143 11.849790  0.0000
## minutter         0.204885 0.2698103 143  0.759366  0.4489
## minutter:treatA -1.515154 0.3882902 143 -3.902118  0.0001
##  Correlation: 
##                 (Intr) minttr
## minutter        -0.201       
## minutter:treatA  0.002 -0.667
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -1.32425075 -0.63786627 -0.04937682  0.64238165  1.97106796 
## 
## Number of Observations: 161
## Number of Groups: 16

ggplot(compute_residuals(model6a), aes(sample = large_resids_scaled)) +
    geom_qq() +
    geom_qq_line()

ggplot(compute_residuals(model6a), aes(x = fitted_values, y= large_resids_scaled)) +
    geom_point() + 
    geom_smooth()

ggplot(compute_residuals(model6a), aes(sample = conditional_resids_scaled)) +
    geom_qq() +
    geom_qq_line()

plot_df <- model6a$data %>% 
    mutate(pred_fix = predict(model6a, level = 0), # fixed effects
           pred = fitted(model6a)) # patient-level effects

ggplot(plot_df, aes(x = minutter, y = pred_fix, colour = treat)) +
    geom_line() +
    geom_point()

ggplot(plot_df, aes(x = minutter, y = pred, colour = patient)) +
    geom_line() +
    geom_point() +
    facet_wrap(~ treat, nrow = 1)

Compare with model 5A

model6a_ml = update(model6a, method = "ML")

anova(model5a_ml, model6a_ml)

##            Model df      AIC      BIC    logLik   Test L.Ratio p-value
## model5a_ml     1 24 620.8233 694.7770 -286.4117                       
## model6a_ml     2  6 602.6345 621.1229 -295.3173 1 vs 2 17.8112  0.4682

# The estimated effect at 5 minutes (because linear we multiply by 5)
intervals(model6a, which="fixed")$fixed[3,] * 5

##      lower       est.      upper 
## -11.413423  -7.575772  -3.738121

result <- result %>% 
    rbind(c("6a: SP(POW)-model med lineær effekt af tid", "-7.58 (-11.41, -3.74)", 0.0001))

b) Random-effect linear regression

model6b <- lme(pressure ~ minutter + treat:minutter,  
               random = ~ 1 + minutter | patient, 
               na.action = na.omit, 
               data = lang)

summary(model6b)

## Linear mixed-effects model fit by REML
##  Data: lang 
##        AIC      BIC    logLik
##   631.9417 653.3798 -308.9708
## 
## Random effects:
##  Formula: ~1 + minutter | patient
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev    Corr  
## (Intercept) 7.4684951 (Intr)
## minutter    0.7919332 -0.289
## Residual    1.0937369       
## 
## Fixed effects: pressure ~ minutter + treat:minutter 
##                     Value Std.Error  DF   t-value p-value
## (Intercept)     20.761744 1.8739862 143 11.078920  0.0000
## minutter         0.172533 0.2719200 143  0.634499  0.5268
## minutter:treatA -1.564592 0.3981502 143 -3.929652  0.0001
##  Correlation: 
##                 (Intr) minttr
## minutter        -0.226       
## minutter:treatA  0.002 -0.648
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.02054204 -0.48388203 -0.03573534  0.49987358  3.00947584 
## 
## Number of Observations: 161
## Number of Groups: 16

ggplot(compute_residuals(model6b), aes(sample = large_resids_scaled)) +
    geom_qq() +
    geom_qq_line()

ggplot(compute_residuals(model6b), aes(x = fitted_values, y= large_resids_scaled)) +
    geom_point() + 
    geom_smooth()

ggplot(compute_residuals(model6b), aes(sample = conditional_resids_scaled)) +
    geom_qq() +
    geom_qq_line()

ggplot(compute_residuals(model6b), aes(x = minutter, y = conditional_resids_scaled, colour = patient)) +
    geom_line() +
    geom_point() +
    facet_wrap(~ treat, nrow = 1)

plot_df <- model6b$data %>% 
    mutate(pred_fix = predict(model6b, level = 0), # fixed effects
           pred = fitted(model6b)) # patient-level effects

ggplot(plot_df, aes(x = minutter, y = pred_fix, colour = treat)) +
    geom_line() +
    geom_point()

ggplot(plot_df, aes(x = minutter, y = pred, colour = patient)) +
    geom_line() +
    geom_point() +
    facet_wrap(~ treat, nrow = 1)

# Again, linear model so we just multiply the coefficient by 5 
intervals(model6b, which = "fixed")$fixed[3, ] * 5 

##      lower       est.      upper 
## -11.758061  -7.822959  -3.887857

result <- result %>% 
    rbind(c("6b: Random-effects regression", "-7.82 (-11.76, -3.89)", 0.0001))

Models 6a and 6b have the same mean-value structure, so we can compare them directly. The AR(1) model has a lower AIC so we prefer that as it fits the data better.

anova(model6a, model6b)

##         Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## model6a     1  6 601.4561 619.8316 -294.7280                        
## model6b     2  7 631.9417 653.3798 -308.9708 1 vs 2 28.48562  <.0001

Q7

Including the outlier doesn’t change the overall conclusion (active-arm patients have a steeper trend) and the trend is even amplified. This makes sense if we consider the data of that outlier.

model7 = lme(pressure ~ minutter + treat:minutter, 
             random = ~ 1 | patient, data = lang5,
             correlation = corCAR1(0.9, ~ minutter | patient, FALSE), 
             na.action = na.omit)

summary(model7)

## Linear mixed-effects model fit by REML
##  Data: lang5 
##        AIC      BIC    logLik
##   661.4622 680.2059 -324.7311
## 
## Random effects:
##  Formula: ~1 | patient
##         (Intercept) Residual
## StdDev:    5.057275  4.41461
## 
## Correlation Structure: Continuous AR(1)
##  Formula: ~minutter | patient 
##  Parameter estimate(s):
##       Phi 
## 0.9011342 
## Fixed effects: pressure ~ minutter + treat:minutter 
##                     Value Std.Error  DF   t-value p-value
## (Intercept)     20.362616 1.6271077 152 12.514608  0.0000
## minutter         0.192781 0.2812642 152  0.685409  0.4941
## minutter:treatA -1.614489 0.3907038 152 -4.132258  0.0001
##  Correlation: 
##                 (Intr) minttr
## minutter        -0.209       
## minutter:treatA  0.002 -0.689
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -1.42325828 -0.57640556  0.02810233  0.55744791  1.91974344 
## 
## Number of Observations: 171
## Number of Groups: 17

intervals(model7)$fixed[3, ] * 5

##      lower       est.      upper 
## -11.932000  -8.072444  -4.212889

result <- result %>% 
    rbind(c("7: som 6a, med person nr. 13", "-8.07 (-11.93, -4.21)", 0.0001))

result

## # A tibble: 8 x 3
##   metode                                     behandlingseffekt_5_minutt… p_value
##   <chr>                                      <chr>                       <chr>  
## 1 3d: T-test efter 5 minutter                -0.27 (-7.57, 7.01)         0.92   
## 2 3e: T-test på ændring 0-5 minutter         -7.61 (-12.88, -2.33)       0.015  
## 3 3e: ANCOVA, baseline som kovariat          -5.20 (-12.76, 2.37)        0.14   
## 4 4b: CS-model med minutter som class        -8.42 (-11.33, -5.53)       0      
## 5 5: SP(POW)-model med minutter som class    -7.80 (-12.00, -3.61)       3e-04  
## 6 6a: SP(POW)-model med lineær effekt af tid -7.58 (-11.41, -3.74)       1e-04  
## 7 6b: Random-effects regression              -7.82 (-11.76, -3.89)       1e-04  
## 8 7: som 6a, med person nr. 13               -8.07 (-11.93, -4.21)       1e-04

# diagnostics as for model 6a